#include <iostream>
using namespace std;

#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4ThreeVector.hh"
#include "Randomize.hh"

#include "RICHMessenger.hh"
#include "RICHPrimaryGeneratorAction.hh"

RICHPrimaryGeneratorAction::RICHPrimaryGeneratorAction() {
	rayGun = new G4ParticleGun(1);
	G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
	G4ParticleDefinition *particle = particleTable->FindParticle("mu-");
	rayGun->SetParticleDefinition(particle);
	rayGun->SetParticleTime(0);
	rayGun->SetParticlePosition(G4ThreeVector(0, 0, 1.5 * m));
	rayGun->SetParticleMomentumDirection(G4ThreeVector(0, 0, -1));
	rayGun->SetParticleEnergy(.7 * GeV);
}

RICHPrimaryGeneratorAction::~RICHPrimaryGeneratorAction() {
	delete rayGun;
}

void RICHPrimaryGeneratorAction::GeneratePrimaries(G4Event *event) {
	/*
	 rayGun->SetParticlePosition(
	 G4ThreeVector(-9.9 * cm + 19.8 * cm * CLHEP::RandFlat::shoot(),
	 -9.9 * cm + 19.8 * cm * CLHEP::RandFlat::shoot(),
	 -1.49 * m));
	 */
	rayGun->GeneratePrimaryVertex(event);
}

double RICHPrimaryGeneratorAction::getCosmicMomentum(double randFlat) {
	double a = 0.00559394, b = 0.54543, c = 0.561472, d = 0.00936388, D =
			-0.402605, z = 0.140012, pi = 3.14159265;

	double arg = d / a * (randFlat + D);
	double gamma = log(1 - arg * arg);
	double alpha = 2 / pi / z + gamma / 2;

	return exp(
			sin(arg) / b * sqrt(sqrt(alpha * alpha - gamma / z) - alpha)
					- c / b) * GeV;
}
